The frustrated Heisenberg antiferromagnet on the honeycomb lattice: 
A candidate for deconfined quantum criticality 
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We study the ground-state (gs) phase diagram of the frustrated spin-i J1-J2-J3 antiferromagnet 
with J2 = Jz = kJi on the honeycomb lattice, using coupled-cluster theory and exact diagonaliza- 
tion methods. We present results for the gs energy, magnetic order parameter, spin-spin correlation 
function, and plaquette valence-bond crystal (PVBC) susceptibility. We find a Neel antiferromag- 
netic (AFM) phase for k, < k c1 ~ 0.47, a collinear striped AFM phase for k, > /t c , ~ 0.60, and a 



paramagnetic PVBC phase for k c 



The transition at k C2 appears to be of first-order 



type, while that at k c1 is continuous. Since the Neel and PVBC phases break different symmetries 
our results favor the deconfinement scenario for the transition at k c , . 



PACS numbers: 75.10.Jm, 75.50.Ee, 03.65.Ca 



Frustrated quantum Heisenberg magnets are 
paradigms of those systems that may be used to 
study quantum phase transitions (QPTs) between qua- 
siclassical magnetically ordered phases and magnetically 
disordered quantum phases. In particular, if the quasi- 
classical phase and the quantum phase spontaneously 
break different symmetries, the Landau-Ginzburg- 
Wilson scenario of continuous phase transitions does 
not hold and the concept of deconfined quantum crit- 
icality [l[ becomes relevant. A canonical model for 
studying such QPTs is the spin-i Heisenberg anti- 
ferromagnet (HAFM) with antiferromagnetic (AFM) 
nearest-neighbor (nn) Ji coupling and frustrating AFM 
next-nearest-neighbor (nnn) J2 coupling (namely the 
J1-J2 model) on the square lattice (see, e.g., Refs. 
0-0] )■ It is now widely acknowledged that this model 
undergoes a transition from a quasiclassical Neel state 
at J2/J1 < 0.4 to a quantum paramagnetic (QP) phase, 
which is not magnetically ordered, for 0.4 < J2/J1 <■ 0.6, 
and thence to a collinear striped phase for J2/J1 > 0.6. 
The synthesis of layered magnetic materials that might 
be described by the square-lattice J±- Ji model [1, Q has 
also encouraged further theoretical studies of this model. 
Although the square- lattice spin-i J1-J2 HAFM has 
been intensively studied for over 20 years, no consensus 
yet exists on such fundamental issues as the nature of 
the QP phase and of the type of transition into it. In 
particular, there is still controversy over whether the 
scenario of deconfined criticality holds for this model 
(see, e.g., Refs. [HH). 

Related magnetic systems may shed light onto this con- 
troversy and one strongly related model that has received 
much attention recently has been the HAFM on the hon- 
eycomb lattice. One reason for this interest is that a 
spin-liquid phase has been found for the exactly solvable 
Kitaev model in which the spin-i particles reside on 



a honeycomb lattice. The honeycomb lattice is also rele- 
vant to the very active research field of graphene, where 
Hubbard-like models on this lattice may be appropriate 
to describe the relevant physics [llj . Interestingly, Meng 
et al. [l2[ have very recently found clear evidence that 
the quantum fluctuations are strong enough to establish 
an insulating spin-liquid phase between the non-magnetic 
metallic phase and the AFM Mott insulator for the Hub- 
bard model on the honeycomb lattice at moderate values 
of the Coulomb repulsion U . At very large values of U, 
the latter phase corresponds to the HAFM on the bi- 
partite honeycomb lattice that contains a ground-state 
showing Neel long-range order (LRO). However, higher 
order terms in the t/U expansion of the Hubbard model 
may lead to frustrating exchange couplings in the cor- 
responding spin-lattice model where the HAFM with nn 
exchange is the leading term in the large- U expansion. 
This unexpected result [12| has also stimulated interest 
in the frustrated HAFM on the honeycomb lattice (see 
e.g. Refs. [H-[l6j]). Indeed, the available literature sug- 
gests a frustration-induced QP phase for the frustrated 
spin-i HAFM on the honeycomb lattice 13-lj|. These 
theoretical findings are consistent with the recent exper- 
imental observations that the spin-| honeycomb-lattice 
HAFM Bi 3 Mn 4 0i2(N0 3 ) shows a spin-liquid-like behav- 
ior at temperatures much lower than the Curie- Weiss 
temperature [20( | . 

One approach to shed additional light on the QP re- 
gions of such J\-Ji models has been to extend its param- 
eter space by the inclusion of next-next-nearest-neighbor 
(nnnn) couplings J3 as well. This yields the J1-J2-J3 
model (see, e.g., Ref. [2l[ and references cited therein), 

H = Ji Si ■ Sj + J 2 ^2 Si-Sfc + J3 ^2 s *' Sfe ' w 

(id) ((i,k)) «<U)» 

where i runs over all lattice sites on the lattice, and j 



2 



runs over all nn sites, k over all nnn sites, and I over 
all nnnn sites to i, respectively, counting each bond once 
and once only. Each site i of the lattice carries a spin-i 
particle with spin operator Sj. 

Here we wish to study the spin-i J\-Ji~Jz HAFM 



model of Eq. ([T]) on the honeycomb lattice 13l4l9l| . The 
lattice and the exchange bonds are illustrated in Fig. Q] 
(right panel), where the quasiclassical Neel and collinear 
striped phases are also shown. Classically, this system ex- 
hibits the two collinear phases illustrated in Fig. [1] (right 
panel) as well as a spiral phase. These phases meet in a 
triple point at J2 = J3 = Ji/2, (see, e.g., Refs. fl3l flit). 
Henceforth we set Jj = 1 and we further restrict our- 
selves to the case J2 = J3 = kJ\. The motivation to fo- 
cus on J2 = J3 also comes from the classical case, where 
the honeycomb-lattice model and the square-lattice J\- 
J2 HAFM model have similar ground-state properties, 
namely: (i) there is a direct transition between the Neel 
and the collinear striped states at k = 0.5; and, (ii) the 
classical ground state is highly degenerate at the transi- 
tion point. Hence, the two systems are very similar clas- 
sically, and so the honeycomb system is a likely candidate 
to exhibit a QP phase and possibly also a deconfined crit- 
ical point between this QP phase and the quasiclassical 
phases. Moreover, we note again that non-zero J3 bonds 
are likely from higher-order terms in the t/U expansion 
of the Hubbard model at moderate U values [l2| . 

The treatment of highly frustrated quantum magnets 
is notoriously challenging, since only a relatively few ac- 
curate methods exist. New promising approaches such as 
projected entangled pair states are not yet accurate 
enough. The finite-size extrapolation of numerical exact 
data for finite-lattice systems, used successfully for the 
square-lattice model @, 13| , is also somewhat less efficient 
for the honeycomb lattice, since (i) the unit cell con- 
tains two sites, (ii) nnnn couplings J3 are included, and 
(iii) there are only a few finite lattices with full point 
group symmetry [18|. Naturally, quantum Monte Carlo 
methods are severely limited in the presence of frustra- 
tion due to the infamous "sign problem." The coupled- 
cluster method (CCM) when evaluated to high orders of 
approximation provides a reasonably accurate approach 
to determine the position of QPT points 
as well as providing evidence on the nature of such QP 
phases [|. Here we use this method together with com- 
plementary results provided by the exact diagonalization 
(ED) technique for a finite lattice of N = 32 sites. 

The CCM is an intrinsically size-extensive method that 
provides results in the limit N — ¥ 00 from the outset. 
The CCM requires us to input a model (or reference) 
state [27T - I29] . with respect to which the quantum corre- 
lations are expressed. Here we use the Neel and striped 
model states shown in Fig. Q] (right panel). We use also 
the very well-tested localized LSUBm truncation scheme 
in which all multi-spin correlations are retained in the 
CCM correlation operators over all distinct locales on 
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1: (Color online) Left: Results for the gs energy, 
E/N (Ji = 1 and J 3 = J 2 ) via ED (N = 32), CCM 
LSUBm (m = 6, 8, 10, 12) and extrapolated m — > 00 
(see text). Right: The J1-J2-J3 honeycomb model; for 
the Neel and the collinear striped model states. The 
arrows represent spins located on specific lattice sites 
[symbolized by •]. 



the lattice defined by m or fewer contiguous sites. The 
method of solving for higher orders of LSUBm approx- 



27, m 



The 



imations is discussed in detail in Rcfs. 
number of independent spin configurations taken into ac- 
count in the correlation operator increases rapidly with 
the truncation index m. For example, at the highest 
level of approximation used here, namely LSUB12, we 
take into account 750490 such configurations. To ob- 
tain results in the exact m — > 00 limit the raw LSUBm 
data must be extrapolated. Although there are no ex- 
act extrapolation rules, a great deal of experience has 
been garnered by now for the gs energy and for the 
magnetic order parameter (sublattice magnetization) M. 
For the gs energy per spin a well-tested and very ac- 
curate extrapolation ansatz (see, e.g., Refs., (0, [231 — 



2a|2J,|30j) is E(m)/N = ao + a^- 2 +a 2 m,- 4 , while for 
the magnetic order parameter M different schemes have 
been employed for different situations 3l(. Here we use 
M(m) = bo + bimT 1 / 2 + &2«i ~ 3 ^ 2 - This scheme has been 
found to be appropriate for systems showing a gs order- 
disorder QPT (see, e.g., Ref. 0, ED, H S|). Since the 
hexagon is an important structural element of the hon- 
eycomb lattice we decided to take into account for the 
extrapolations only LSUBm data with m > 6. 

CCM results for the gs energy per spin E/N are shown 
in Fig. []] using both the Neel and collinear striped model 
states. They are clearly well-converged for all values of J2 
shown. The corresponding extrapolated LSUBoo results 
are also shown in the regimes where real solutions exist 
for the entire data set, for each choice of reference state. 
These results compare well to those from ED calculations 
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FIG. 2: (Color online) CCM results for the gs order 
parameter, M ( Ji = 1 and J3 = J2) for the Neel and 
striped phases using LSUBm (m = 6, 8, 10, 12) and 
extrapolated m —> 00 (see text). 

for N = 32, also shown in Fig. [U The CCM results 
show a clear intermediate region around k ~ 0.5 in which 
neither of the quasiclassical AFM states is stable. The 
ED results demonstrate a "kink- like" behavior in E/N 
at k « 0.6, (cf. Refs. @,[H), which is a first hint of a 
first-order QPT (and see the discussion below). 

Wc now discuss the magnetic order parameter M in or- 
der to investigate the stability of quasiclassical magnetic 
LRO. CCM results for M are shown in Fig. [U The ex- 
trapolated Neel order parameter vanishes at k Ci 0.466, 
whereas the corresponding collinear striped order param- 
eter goes to zero at k C2 ss 0.601. These values can 
be considered as CCM estimates for the QCPs of the 
model [32[. They arc in good agreement with those re- 
cently obtained by the Schwinger-boson approach [liij . 
the functional renormalization group approach (l6| . and 
by an ED approach [Hj]. From the shape of the order- 
parameter curve in Fig.[2]it seems probable that the tran- 
sition at k Ci is continuous, whereas the very steep fall 
near k C2 may indicate a first-order transition. This ob- 
servation is supported by the spin-spin correlations func- 
tions shown in Fig. [3l which also show a similar "smooth" 
versus "steep" behavior near the points k Ci and k C2 , re- 
spectively. Moreover, the CCM and ED data presented in 
Fig.[3]arc in good quantitative agreement, which strongly 
supports the conclusions that have been drawn above re- 
garding the likely nature of the QPTs. 

Hitherto we have gathered strong evidence about the 
region of the intermediate nonmagnetic quantum phase 
by using the high-order CCM and ED techniques. How- 
ever, the question remains as to what is the nature 
of this phase. We note that numerical evidence was 
found very recently of a plaquette valence-bond crys- 
tal (PVBC) phase near k = 0.5 jljj. Note that such 
a PVBC phase is a strong candidate for the QP phase 
of the J 1 J 2 model on the square lattice . In or- 
der to analyze the possibility of such a PBVC phase 
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FIG. 3: (Color online) CCM LSUB12 and N = 32 ED 
results for the spin-spin correlation function (s(0) • s(r)) 
( Ji = 1 and J3 = J2). Values r = 1, 1.732, 2 correspond 
respectively to nn, nnn, and nnnn pairs of spins. 



we first consider a generalized susceptibility xf that de- 
scribes the response of the system to a "field" opera- 
tor F (see, e.g., Refs. 0, HI). We thus add a field term 
F = 5 O to the Hamiltonian ([T]), where O is an opera- 
tor that in our case corresponds to the possible PVBC 
order illustrated in Fig. [4] (right panel), and which hence 
breaks the translational symmetry of H . Wc now cal- 
culate the energy per site E(8)/N = e(S) for the per- 
turbed Hamiltonian H + F, by using the CCM for both 
the Neel and the collinear striped reference states, and 
define the susceptibility as \f = — (d 2 e(S))/ (dS 2 )\ s=Q . 
Rather less empirical experience is available regarding 
the extrapolation of the CCM LSUBm data for \f than 
for other quantities such as the gs energy or M. How- 
ever, we have tested several extrapolation schemes and 
found that the extrapolation used for the gs energy, i.e. 
Xi?(m) = Co + cito~ 2 + C2TO -4 , fits the data most ac- 
curately. For a crosscheck of this extrapolation scheme 
we also use corresponding extrapolation of the inverse 
susceptibility, i.e. X^ 1 (m) = do + d^rnr 2 + d2in~ 4 , and 
find consistent results as shown in Fig. 2] An instabil- 
ity of the gs against the perturbation F is indicated by 
a divergence of \f or equivalently a zero point of Xf 1 ■ 
This fact might favor slightly the extrapolation of x^ 1 
when searching for such instabilities. Our results are in 
accordance with those of Ref. [l5| and they clearly fa- 
vor a PVBC intermediate quantum phase instead of a 
structureless spin-liquid phase. The extrapolated inverse 
susceptibility vanishes on the Neel side at k w 0.473 and 
on the collinear striped side at k w 0.586. These val- 
ues are in very good agreement with the critical values 
k Ci ks 0.466 and k C2 ps 0.601 already found by the CCM 
order parameter. Thus, it is most likely that the PVBC 
phase occurs at (or is very close to) the transition points 
where the quasiclassical magnetic LRO breaks down. Al- 
though there is a steep downturn of Xf 1 a t K c 2 > wc find a 
smooth behavior for Xf 1 & t K c ± ■ This evidence strongly 
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FIG. 4: (Color online) Left: Results for the inverse 
plaqucttc susceptibility, 1/x (Ji = 1 and J3 = J 2) via 

CCM LSUBm (m = 6, 8, 10, 12), and extrapolated 
(m — > 00) for both x an d 1/x ( see text). Right: The 
perturbations (fields) F = 5 O for the plaquette 
susceptibility x- Thick red and thin black lines 
correspond respectively to strengthened and weakened 
nn exchange couplings, where O = a ij s i ' s ji an( i 

the sum runs over all nn bonds, with a%j = +1 (—1) for 
thick red (thin black) lines. 



supports our conclusion drawn from Figs. [TJ [21 and [3] that 
the phase transition is of first-order type at k C2 while that 
at K C1 is of continuous type. 

The possibility of deconfined quantum criticality for 
the frustrated honeycomb HAFM was first pointed out 
in Ref. [l[. The first evidence supporting such a sce- 
nario for the J1-J2-J3 model on the honeycomb lattice 
has been found very recently in Ref. [15j . Our CCM re- 
sults are seen to be highly converged and to agree well 
with those results of the ED method. Our results for the 
gs energy, magnetic order parameter, spin-spin correla- 
tion function, and plaqucttc susceptibility all point to 
collinear phases separated by a magnetically disordered 
phase for 0.47 < k < 0.60. We have shown that the na- 
ture of this intermediate phase is likely to be of plaqucttc 
valence-bond crystalline (PVBC) type. Although we find 
indications of a first-order transition at k c , 2 rj 0.60, the 
transition between the Neel and the PVBC phases at k Ci 
appears continuous. Since the Neel and the PVBC phases 
break different symmetries and we have shown that they 
are likely to meet at k c1 ~ 0.47, our results present strong 
evidence that the frustrated spin-i HAFM on the hon- 
eycomb lattice contains a continuous QPT described by 
the scenario of deconfined quantum criticality. 
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